library(optparse)


# 1 -----------------------------------------------------------------------

option_list <- list(  #构建参数列表
  make_option( c("-d", "--data"), type = "character", default = FALSE, help = "输入文件"),  
  make_option( c("-c", "--cutoff"), type = "integer", default = 0.05, help = "选取阈值"),
  # type为输入参数的类型，有“logical”, “integer”, “double”,“complex”, “character”这几种类型
  # default为该参数的默认值，当用户没有输入该值时使用默认值，
  #help为输出的帮助信息
) #有几个参数就用几个make_option
opt_parser = OptionParser(option_list=option_list)  #解析参数
opt = parse_args(opt_parser)   #解析参数 
x1 = opt$data   #将输入的-d参数传递给x
y1 = opt$cutoff  #将输入的-c参数传递给y

Afunction<- function(x,y){
  
}   #函数本身

write.table(Afunction(x1,y1),"result.csv",quote = F,col.names = T,row.names = F,sep=',')  #调用函数并将结果保存成文件输出

# DESeq2数据准备 --------------------------------------------------------------
# 需要两个文件
# 矩阵文件 aggr_
#
DESEq2_gene <- aggr_TCGA_count
y <- rownames(aggr_TCGA_count)
type <- unlist(lapply(y,function(y) strsplit(as.character(y),"-")[[1]][4]))
aggr_TCGA_count$type <- type
table(aggr_TCGA_count$type)

aggr_TCGA_count <- na.omit(aggr_TCGA_count)
aggr_TCGA_count$id <- rownames(aggr_TCGA_count)
sum_count <- merge(aggr_TCGA_count,LUNG_survival,by = 'id')
for(i in 1:58389){
  sum_count[,i] <- as.numeric(sum_count[,i])
}
LUNG_survival <- sum_count[,58388:58389]
y <- rownames(LUNG_survival)
type <- unlist(lapply(y,function(y) strsplit(as.character(y),"-")[[1]][4]))
LUNG_survival$type <- type

group <- factor(LUNG_phenotype_01A$M,levels = c("M0","M1"))

DESEq2_design<-model.matrix(~-1+group)

rownames(DESEq2_design) <- rownames(LUNG_phenotype_01A)
aggr_TCGA_count <- sum_count[,1:58387]
DESEq2_gene <- aggr_TCGA_count
DESEq2_gene <- 2^DESEq2_gene-1
DESEq2_gene <- as.data.frame(t(DESEq2_gene))
DESEq2_gene <- round(DESEq2_gene,0)

# DESeq2 ------------------------------------------------------------------
# 表达矩阵 DESEq2_gene，设计矩阵 DESEq2_design
library('DESeq2')#加载包
all(rownames(DESEq2_design)==colnames(DESEq2_gene))#确保表达矩阵的列名与分组矩阵行名相一致

colData <- data.frame(row.names = colnames(DESEq2_gene), group)


dds <- DESeqDataSetFromMatrix(countData=DESEq2_gene, colData=colData, design<- ~ group)  #DESeq2数据格式的构建
#如果做过log转换，需要还原为counts
#express_rec <- 2^express_rec-1
head(dds)
dds <- dds[ rowSums(counts(dds)) > 1, ] #过滤一些low count的数据；
dds <- DESeq(dds)#DESeq进行标准化；
resultsNames(dds)
res <- results(dds)
summary(res)#查看经过标准化矩阵的基本情况；
mcols(res,use.names = TRUE)
res_data <- merge(as.data.frame(res),as.data.frame(counts(dds,normalize = TRUE)),by = 'row.names',sort = FALSE)
res_data <- res_data[,1:7] 
dif_DESeq2<-res_data[res_data[,"padj"]<0.05,]#pvalue or padj
dif_DESeq2 <- dif_DESeq2[abs(dif_DESeq2[,"log2FoldChange"])>1.2,]
dif_DESeq2 <- na.omit(dif_DESeq2)

write.csv(dif_DESeq2,file = "result/DESeq2/dif.csv",quote = F)
write.csv(res_data,file = "result/DESeq2/dif_all.csv",quote = F)